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Finite Difference Time Marching 
in the Frequency Domain: 

A Parabolic Formulation for the 
Convective Wave Equation 

An explicit finite difference iteration scheme is developed to study harmonic sound 
propagation in ducts. To reduce storage requirements for large 3D problems, the 
time dependent potential form of the acoustic wave equation is used. To insure that 
the finite difference scheme is both explicit and stable, time is introduced into the 
Fourier transformed ( steady-state ) acoustic potential field as a parameter. Under a 
suitable transformation, the time dependent governing equation in frequency space 
is simplified to yield a parabolic partial differential equation, which is then marched 
through time to attain the steady-state solution. The input to the system is the ampli- 
tude of an incident harmonic sound source entering a quiescent duct at the input 
boundary, with standard impedance boundary conditions on the duct walls and duct 
exit. The introduction of the time parameter eliminates the large matrix storage 
requirements normally associated with frequency domain solutions, and time 
marching attains the steady-state quickly enough to make the method favorable when 
compared to frequency domain methods. For validation, this transient-frequency 
domain method is applied to sound propagation in a 2D hard wall duct with plug 
flow. 


Introduction 

Both steady-state (frequency domain) and transient (time 
domain) finite difference and finite element techniques have 
been developed to study sound propagation in ducts. To date, 
the numerical solutions have generally been limited to moderate 
frequency sound and mean flow Mach numbers. Wavelength 
resolution problems have prevented a broader range of applica- 
tions of the numerical methods. A fine grid is required to resolve 
the short wavelengths associated with high frequency sound 
propagation with high Mach numbers. 

The present paper focuses on sound propagation in a duct 
with flow. This is the first step in a larger research effort to 
develop efficient numerical techniques to predict high frequency 
sound propagation around 3D aircraft inlet nacelles with large 
subsonic inlet Mach numbers. The paper begins with a descrip- 
tion of the problem, an explanation of why the transient ap- 
proach is employed, and a description of the governing equa- 
tions. The bulk of the paper describes the development of a 
stable, explicit finite difference scheme by a transformation of 
the governing hyperbolic wave equation. The scheme is iterated 
in time to converge to the steady-state solution associated with 
a Fourier transform solution. 

Problem Statement 

The problem under consideration here is the steady-state 
propagation of sound, represented by the perturbation acoustic 
potential, through a 2D rectangular duct with flow. The goal of 
the paper is to develop a stable, explicit finite difference scheme 
that incorporates an impedance condition applied at the duct 
exit and rigid body boundary conditions on the duct walls. For 
the numerical examples considered here, the grid system shown 
in Fig. 1 is employed. 


Contributed by the Technical Committee on Vibration and Sound for publica- 
tion in the Journal of Vibration and Acoustics. Manuscript received Feb. 
1995; revised July 1995. Associate Technical Editor: G. H. Koopmann. 


Transient Approach 

In frequency domain approaches, pressure and acoustic ve- 
locity are assumed to be harmonic functions of time. The matri- 
ces associated with numerical solutions of the governing equa- 
tions in the frequency domain have extremely large storage 
requirements. In similar 3D electromagnetic applications, fre- 
quency domain approaches are limited to several hundred thou- 
sand unknowns (still a small problem in 3D). Larger problems 
encounter roundoff errors and conditioning problems that pre- 
vent a reliable solution. 

On the other hand, an explicit transient method generates no 
matrices, and is generally faster than a frequency domain ap- 
proach (Baumeister, 1980a). Miller has developed explicit rela- 
tionships showing the time advantage of transient over steady- 
state techniques (Miller, 1988). Consequently, the method of 
choice in the present paper is a time dependent method. 

Governing Equations 

The governing differential equations for studying acoustic 
propagation in ducts can be formulated in terms of the constitu- 
tive continuity and momentum equations or in terms of potential 
flow. The constitutive equation approach can handle a general 
3D sheared flow while the potential approach is limited to 3D 
inviscid flow. The major advantage of the potential flow ap- 
proach over the constitutive equation approach comes in de- 
creased storage requirements associated with only one depen- 
dent variable. 

Fortunately, acoustic propagation can often be reasonably 
modeled by an inviscid approximation. For single mode JT15D 
engine inlet nacelles data, a previous finite element study 
(Baumeister and Horowitz, 1984) employing the potential for- 
mulation in the frequency domain showed good agreement with 
experimental data — in the far field radiation pattern as well as 
suppressor attenuation. Due to this success, the problem under 
consideration here is formulated in terms of an acoustic poten- 
tial. 


622 / Vol. 118, OCTOBER 1996 


Transactions of the ASME 



Fig. 1 Structured FD-TD mesh for rectangular duct 


For inviscid, nonheat conducting and irrotational flow, the 
linearized potential wave equation can be written as 

0 =/ 2 <K - (c 2 - § 2 )<^ - (c 2 - $ 2 )<^ 


and 


c 2 = 1 - 3(7 - l)Mj ( 6 ) 


ip' = -/</>; -MfM ( 7 ) 

p 

Equations ( 1 ) and (5) are the basic equations used to estab- 
lish a finite difference formulation. However, care must be taken 
when discretizing derivative terms to insure that the resulting 
scheme is both stable and explicit. Note that it is easy to develop 
a stable implicit method, but this would yield a matrix formula- 
tion that is no better than a frequency domain approach. It 
turns out that the proper treatment of the mixed time and space 
derivative terms [which appear on the second line of Eq. ( 1 )] 
is critical for maintaining stability. 

The implicit formulation is obtained by approximating the 
mixed partials by 


J. r 1 1 ±tk • ±fk , jl r fc— 1 JL / fc+ 1 A\ r 1 

, _ <Pi.j + + <f>i+u + <Pi,j - 2<pij - <Pj-u - <Pi+i.j 




2 Ax At 
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+ 2 ( 5 ^ + + 2 ($ x % + 

-(7-1 )(/0; + ^i + $,0; )($„ + $„,) (i) 


( 8 ) 
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1 + <Kf 1 + + 4>p_ x - 2(f> lj - &jtt - <f>'J: 


2AyAt 


(9) 


with 

c= {1 -|( 7 - 1)($ 2 + $ 2 )} 1/2 


where i and j denote the space indices for the nodal system 
shown in Fig. 1, fc is the time index defined by 

t k+l =t k + At (10) 


where <f>(x, 7 ) represents the steady mean flow potential and 
< p'(x , 7 , t) the acoustic potential. The relationship between 
acoustic pressure P' and mean and acoustic potential is 

z p; = -f<i>' t '-($si> , x + $ y <f>' y )' o) 

p 

For the special case of plug flow, the mean flow terms in Eq. 
( 1 ) become 

% = % = = 0 § x = M f (4) 

Substituting Eq. (4) into Eqs. (1), (2) and (3) yields 

0 = f 2 4> ’ - (c 2 - M]W~ - c 2 <f>'yy + 2 (5) 


and Ax, Ay, and At are the grid spacings. 

Baumeister (1980b) showed that Eq. ( 8 ) can be used in an 
explicit fashion for ID plug flow problems in a duct. However, 
if the flow is not one dimensional or if the region exterior to 
the duct is included, the scheme must be implicit. Fortunately, 
this problem can be circumvented by transforming the potential 
and consequently modifying the governing equation. The details 
follow in the next section. 

Transformation to Frequency Space 

There are several ways to develop a frequency domain formu- 
lation for the general 2D acoustic wave Eq. ( 1 ) or the plug flow 
simplification Eq. (5). The Fourier Transform can be applied if 
the potential has a multi-frequency content. In the monochro- 
matic case, this is equivalent to assuming that 


Nomenclature 

c = steady dimensionless speed of 
sound, C # /Co, Eq. (2) 

D # = dimensional duct height or diame- 
ter 

D = dimensionless duct height, D = 1 
d. = parameter, Eq. (33) 

/ # = dimensional frequency 
/ = dimensionless frequency, 
/ # D # /Co, Eq. (38) 
g = parameter, Eq. (33) 
h = param eter, Eq. (33) 
i = JZ 1 

L = length of duct, L # /D # , Fig. 1 
M f = Mach number at duct entrance 
n — unit outward normal 
P = dimensionless fluid pressure, 


P r = acoustic pressure fluctuation, Eq. 
(3) 

t — dimensionless time, / # f # 
t T = total dimensionless calculation 
time 

At = time step 

Us = dimensionless acoustic velocity, 
Eq. (37) 

x = dimensionless axial coordinate, 
x # /D* 

Ax = axial grid spacing 
7 = dimensionless transverse coordi- 
nate, 7 *!D* 

Ay = transverse grid spacing 

Z # = impedance, Eq. (26) 

7 = ratio of specific heats 
£ = dimensionless impedance, 
Z*/p*oC* 0 , Eq. (26) 


p = steady fluid density, p # / pt 
$ = dimensionless time dependent flow 
potential, $ # /CoD # 

$ = steady mean flow potential, Eq. ( 1 ) 
<f>' = transient acoustic potential, Eq. (1) 
<f> = transient acoustic potential in fre- 
quency space, Eq. (13) 
d = dimensionless frequency, 27 r/ 

Subscripts 

i — axial index 
j = transverse index 
o = ambient or reference condition 

Superscripts 

# = dimensional quantity 
fc = time step 
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4>'(x, y, t ) = i f/(x, y)e~ ,u * r * = 0(x, y)e~ ,2nt (11) 
which, in the case of plug flow [ from Eq. (5)], yields 

0 = (c 2 - M J'jiffxx + c 2 tf/yy + u 2 ifj + (12) 

This equation would be solved numerically using a linear 
system of equations. However, the associated matrix is not posi- 
tive definite, which can lead to numerical difficulties, and which 
precludes the use of iterative techniques (see TRANSIENT 
APPROACH). Therefore, it is desirable to develop an explicit 
finite difference scheme to avoid the use of matrices. In time- 
dependent form, Eq. (1 ) or (5) cannot easily be discretized in 
such a way that the resulting finite difference scheme is both 
stable and explicit in the presence of flow, although it is possible 
to obtain reasonable results in the no-flow case. 

The resolution of these difficulties is achieved by modifying 
the monochromatic transformation (11) to 

y, t) =.4>(x,y, = 0(x, y, t)e~ l27 " (13) 

This differs . frpm the--Glassical monochromatic transformation 
in that the amplitude 0 (no prime) is no longer assumed to be 
independent of time, as in Eq. (11). Employing Eq. (13), the 
time derivatives in Eqs. (1) and (5) can be replaced by the 
following relationships: 

4>\ = [~i2ir<t> + Me- 12 " (14) 

4>'*= [-&*<!>, + Me- 1 *" (15) 

<K = -I (00 = (0«, - 2i27T0 r - (27T) 2 0)e- i2 " (16) 


state solution. The mixed space-time derivative term in Eq. 
(18) prevents an explicit finite difference representation of the 
governing equation, and is therefore dropped. Furthermore, the 
hyperbolic time term in Eq. (18) is also dropped to further 
simplify the governing equation. For plug flow, this produces 
a parabolic “wave” equation of the form: 

—2 ifuj<f> t - (c 2 ~ M 2 )4>xx: + c 2 <j>yy + u 2 (f) + i2u>M f (f) x (20) 

This type of differential manipulation is often associated with 
preconditioning of both time dependent partial differential equa- 
tions (Turkel, Fiterman and Leer, 1993) and time independent 
partial differential equations (Turkel and Amone, 1993) to ac- 
celerate convergence to the steady state solution. Here, though, 
the goal is obtaining a stable, explicit scheme. 

In the absence of mean flow, the parabolic wave equation 
takes on the form 

-2ifu)(f> t = 0^ + 4>yy + u> 2 (f) (21) 

The right-hand side of Eq. (21 ) is identified with the Helmholtz 
equation. In the Fourier transform approach, Bayliss, Goldstein, 
and Turkel (1982) have developed an iterative approach to its 
associated matrix; however, the Gaussian elimination approach 
is still the method of choice. In effect, the time iteration proce- 
dure presented in this paper offers a very simple iterative solu- 
tion to this classic equation. 

Initial and Boundary Conditions 

The duct is assumed to be quiescent at time 0, so that the 
initial condition is 


Under this transformation, the general Eq. ( 1 ) and the plug 
flow Eq. (5) become, respectively 

/ 2 0„ - [ i2f u + (y- 1 )/($* + %)}& 

+ 2 f$ x <p x! + 2 f%4> yt = ( C 2 - $ 2 )0xt + (c 2 - $y)<f>yy 

— 2<&/& y 4> xy + u> 2 <f> + i2u)Q x 4> x + i2U)Qy4>y 

- 2($A* + §y%:)4> x - 2($Ay + %$yy)4>y 

+ (7 ~ 1 )(~ 1(^0 + $ X <P X + 5y0y)($ XC + 5 W ) (17) 

/ 2 0„ - 2ifuj<f> t -I- 2 fM f 4> xt 

= (c 2 — M 2 )0^ 4- c 2 4>yy + u 2 <f) + i2uiMf4> x (18) 

To see the relationship between the two transforms [Eqs. 
(11) and (13)], consider, for simplicity, the case of plug flow. 
The monochromatic transform yields Eq. (12), while the tran- 
sient transformation yields Eq. (18). The only difference is in 
the presence of the time derivative terms on the left-hand side 
of Eq. (18). Physically, the time dependence in 0(x, y, t ) is 
caused by assuming that the duct is quiescent at time 0, and 
that the source is turned on at that instant. 

A series of numerical calculations was performed to investi- 
gate the relationship between 0 and 0 as time progresses. It 
was verified that 0 converges to the steady state 0, so that 

lim 0(x, y, t) = 0(x, y) (19) 

f-K» 

At present, a formal proof of convergence is not available. 

Preconditioning 

At present, transient solutions in frequency space are of no 
interest. Therefore, approximations to Eqs. ( 17) or ( 18) can be 
made, by modifying the time derivative terms, as long as the 
following two factors are taken into account: ( 1 ) the resulting 
scheme must be both stable and explicit, and (2) there must be 
at least one time derivative term in the equation to ensure that 
the scheme can be marched through time to obtain the steady- 


0(x, y, 0) = 0 (22) 

As the equation is iterated in time, the solution builds up to the 
steady-state solution. 

Estimating the source mechanism in a duct is a complicated 
process. For example, duct source conditions have been success- 
fully modeled by a Gutin representation in a wind tunnel (Evers- 
man and Baumeister, 1989), by estimating the measured pres- 
sure modal content for a JT15D Turbofan engine (Baumeister 
and Horowitz, 1984) or most recently by an acoustic analogy 
for ducted fans (Farassat et al., 1994). For the validation pur- 
poses of this paper, however, a simple potential source at the 
duct entrance (x = 0) is assumed: 

0'(O, y, t) = l<r ,2 ’ r ' (23) 

and through Eq. (13) to 

0 ( 0 , y, 0=1 ( 24 ) 

The hard wall condition, on the duct walls, is 


V0* n = 0 


(25) 


where n is the unit outward normal. 

To simulate a nonreflective boundary at the duct exit, the 
difference equation is expressed in terms of an exit impedance 
of the form 




z # 

P 


P' 



dx 


(26) 


where Z # is the acoustic impedance. In ducts, this termination 
impedance for a reflection free termination is very restrictive 
(plane waves only). However, this condition is useful as a first 
order termination in the far field from an open ended duct. 
Note that an advantage of the transient-frequency formulation 
presented here is the capability of using the classical impedance 
concept as developed for the frequency domain. 
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Finite Difference Equations 

The finite difference approximations determine the potential 
at the spatial grid points at time steps t k = kAt. Starting from 
the known initial conditions at t = 0 and the boundary condi- 
tions, the algorithm marches the solution out to later times. 

Away from the duct boundaries, each partial derivative in 
Eq. (20) can be expressed as follows: 


jHI tk - 1 

4> t = 

2 At 

, Ku - Wh + 

*** Ax 2 


<f>yy = 


4 bj J+ i ~ 24lj + 1 


Ay 2 


i 4>Uh - tf-u 
<Px 2Ax 

4> - <f>tj 

Substituting these expressions into Eq. (20) yields 


j k+ 

<Pi.j 




At 


2d 2 2c" 


= --=3" 


Ax 2 A 




+ 4>ij + 1 


Ay 2 




(27) 

(28) 

(29) 

(30) 

(31) 


ished in magnitude. This is a desirable property, since the nu- 
merical formulation cannot distinguish between an error and a 
small acoustic mode. 

The von Neumann stability analysis does not take into ac- 
count boundary conditions. For stability, gradient boundary con- 
ditions require the use of smaller At than predicted by Eq. ( 34) . 

Numerical Examples 

The following problem is considered: a plane wave propa- 
gates from the left into a quiescent duct of length one, and the 
acoustic potential field is to be computed in the duct. In the 
three examples that follow, the parabolic transient-frequency 
domain results for plane wave propagation with plug flow are 
compared to the exact results of the steady Fourier transformed 
solutions, given by 

ip(x) = £[<«*/>/<"■*/))* (35) 

For the iterates to converge to ip, a suitable calculation time 
must be specified. Pearson (1953) has shown that the total time 
t T required for the steady-state solution [Eq. (35)] to become 
established in the duct equals the time for a plane wave propa- 
gating at the speed of sound to reach the end of the duct. In 
terms of the dimensional variables, the distance x # a plane wave 
moves in time r # is related to the speed of sound by 


X* = c*t* 

(36) 

so that the dimensionless velocity is 


CJ 

II 

^1- 

. N 

* 

II 

(37) 

where the dimensionless frequency /is 



+ 



+ u 2 4>lj 


(32) 



(38) 


where 

d 2 = c 2 — Mj 
h = iuf 

g = 2iuM f (33) 

Equation (32) is an explicit two step scheme. At t — 0, field 
values at t k ~ 1 are assumed zero because the initial field is quies- 
cent. 

The expressions for the difference equations at the boundaries 
are complicated somewhat by the impedance condition. How- 
ever, a simple integration procedure resolves this problem. 
Baumeister (1980a & 1980b) gives precise details for generat- 
ing the time difference equations at the boundaries. 

Stability 

A von Neumann stability analysis indicates that the method 
is conditionally stable, subject to the condition 



In a typical application, u, /, and M f are set by conditions in 
the duct. Next, the grid spacing parameters Ax and Ay are set 
to accurately resolve the estimated spatial harmonic variation 
of the acoustic field. Finally, At is chosen to satisfy Eq. (34). 

In the von Neumann analysis, conditional stability means that 
the amplification factor, which describes how errors propagate 
from one time step to the next, has magnitude one. Thus, when 
inequality (34) is satisfied, errors are not magnified or dimin- 


According to Pearson, for the steady-state solution ip to de- 
velop and reach the duct exit at x = L, the total calculation 
time must satisfy 

h>- (39) 

Vs 

For the numerical examples below, / and L are each set to 1, 
so t T > 1. 

Pearson also shows that a decaying transient solution propa- 
gates at the speed of sound along with the steady-state solution. 
This transient solution is reflected back into the duct by the exit 
condition. Consequently, some slack time (?r — 3) was added 
in the numerical examples to allow this transient to decay and 
to allow for differences in propagation between the parabolic 
and hyperbolic equations. 

Semi-Infinite Duct. Boundary conditions can introduce im 
stabilities (Baumeister, 1982, andCabelli, 1982) into otherwise 
stable finite difference schemes. Therefore, it is important to 
test the proposed method for convergence in time to the steady- 
state solution in the absence of the exit boundary condition [Eq. 
(26)], and to test independently the effect of the exit boundary 
condition on the solution. Therefore, in this example, the com- 
putational boundary is set at x = 50, far enough away from the 
true boundary x = 1 that any artifacts arising from imperfections 
in the exit boundary condition do not affect the solution in 
[ 0 , 1 ]. 

Two cases are considered — no flow, in Fig. 2(a), and flow 
with Mach number M f = -.5 in Fig. 2(b). For the first case. 
Ax = 0.05, so it would take 1980 time steps ( 1000 forward, 
980 backward) for any artifacts to interfere with the solution. 
Since t T = 5 and At = 0.003, only 1667 steps were run. For 
the second case. Ax = 0.025, so it would take 3960 time steps 
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X axial coordinate 



Fig. 2 Analytical and numerical potential profiles along wall for plane 
wave propagating in a semi-infinite hard wall duct (f = 1). (a) M, = 0 
(Ax = 0.05, At = 0.003, t T = 5.0). (b) M, = -0.5 (Ax = 0.025, At = 
0.001, t T = 5.0). 


(2000 forward, 1960 backward) for any artifacts to interfere 
with the solution. Since t T = 3 and At = 0.001, only 3000 steps 
were run. In both cases, the numerical results show excellent 
agreement with the exact solution. 

Finite Duct, L — 1. In this example, the computational 
boundary is set at the true boundary x = 1 to examine the effect 
of the exit boundary condition [Eq. (26)]. The frequency is 
normalized to 1. Again, two cases were considered — no flow 
(M f - 0) and flow with Mach number M f = -.5. The results 
are shown in Figs. 3(a) and 3(b), respectively. The numerical 
results again match well with the exact results, but it is clear 
that the exit boundary condition does have a slight degrading 
effect on the solution. Notice also that the time step has been 
decreased here, which increases the execution time; however, 
the computational domain is much smaller, which decreases the 
execution time. The total calculation time in this example is t T 
- 4.0. 

Convergence Rate. In this example, the convergence rate 
is studied for the region x = 0 to x = 1 using the semi-infinite 
duct. The results are shown in Fig. 4 for the real and imaginary 


components of the potential and Fig. 5 for the magnitude of the 
potential. As seen in these figures, the numerical solution 
quickly and accurately converges to the exact steady-state solu- 
tion. 

For plane wave propagation, the exact solution to the hyper- 
bolic governing equations indicates the arrival of the “steady 
state” solution at t — 1 [Pearson, 1953, or Baumeister, 1983, 
Eq. (28) renormalized]. As seen in Figs. 4 and 5, the parabolic 
solution does not converge to the steady- state until about time 
t = 2. Thus, the parabolic marching scheme requires more time 
steps to converge than a hyperbolic method. 




Fig. 3 Analytical and numerical potential profiles along wall for plane 
wave propagating in a hard wall duct of unit length and a non-reflecting 
exit condition (f = 1). (a) Mf = 0 (Ax = 0.05, At = 0.0001, t T = 4.0). (b) 
M f = -0.5 (Ax = 0.025, At = 0.00001, t T - 4.0). 
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throughout the domain at the numerical velocity Ax/ At rather 
than the speed of sound. This property may account for the 
slower convergence time for the parabolic formulation as com- 




Fig. 4 Developing history of disturbance propagation in Fourier trans- 
formed domain as a function of time. ( M f = 0, Ax = 0.05, At = 0.003). 
(a) t = 0.25. (b) t = 0.5. (c) t = 0.75. (d) t = 1.00. (e) f = 1.5. (f) f = 2.0. 
(g) t = 2.5. (ft) t = 3.0. 


For parabolic partial differential equations with real coeffi- 
cients, disturbances propagate with infinite speed. In the present 
calculations, numerical solutions to the parabolic Eq. (20) have 
this trait. The finite difference values of the potential propagate 



Fig. 5 Developing history of magnitude of disturbance propagation in 
Fourier transformed domain as a function of time. ( M, = 0, Ax = 0.05, 
Af = 0.003). (a) f = 0.25. (6) t = 0.5. (c) t = 0.75. (d) t - 1.00. (e) t = 1.5. 
(f) t = 2.0. (g) t = 2.5. (h) t = 3.0. 
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Transient solution Transient-frequency Fourier solution 

'Wx.y) solution itfx.y) «Hx,y) 


Fig. 6 Alternate finite difference/element methods in solving wave equations 


pared to the hyperbolic formulation of the appropriate governing 
equations. In the parabolic scheme, acoustic energy is numeri- 
cally transferred ahead of the true wave front. As the algorithm 
proceeds, this transience dies down and the correct steady-state 
solution is attained, using minimal computer storage with run 
times comparable to other methods. 

Historical Perspective 

With the transient-frequency approach developed in this pa- 
per, three different solution techniques are now available to 
solve the wave equation. This is outlined in Fig. 6 for the zero 
mean flow case. 

The Fourier transform of the wave equation was the first 
numerical approach used to study sound propagation in jet en- 
gine ducts (Baumeister and Bittner, 1973). This steady-state 
approach is outlined on the right side of Fig. 6. The governing 
hyperbolic wave equation is transformed to the elliptic Helm- 
holtz equation. Finite difference (FD) and finite element (FE) 
numerical formulations have been employed to solve this equa- 
tion. After the application of the boundary conditions (Fig. 6; 
[BC]), the associated finite difference or finite element global 
matrix is solved for the velocity potential ( or pressure) . Because 
the matrix form of the Helmholtz equation is not positive defi- 
nite, matrix elimination solutions are generally employed. This 
requires extensive storage, as discussed in Transient Approach. 
Conveniently, the steady-state approach allows the direct calcu- 
lation of the potential (or pressure) fields. 


To make the numerical solutions more cost effective, grid 
saving approximations to the governing Helmholtz equation 
have been used (Fig. 6; [Approx.]). Baumeister (1974) em- 
ployed the wave envelope theory. Candel (1986) presents an 
extensive discussion of research accomplishments in this area. 

The transient solution to the wave equation was the second 
numerical approach used to study propagation in ducts, to elimi- 
nate matrix storage requirements, as shown in the left-hand 
column of Fig. 6. Harmonic sound is introduced as a boundary 
condition at the duct entrance into a quiescent duct. FD approxi- 
mations to the wave equation are then solved by an iteration 
process. The calculations are ran until the initial transience dies 
out and steady harmonic oscillations are established. Finally, 
the transient variable 4>' is transformed into the steady-state 
variable if/ associated with the solution of the Helmholtz equa- 
tion. As with the steady-state approach, simplifications to the 
governing equations can reduce computer storage and run times 
(Baumeister, 1986). 

The third option, the transient-frequency technique, is illus- 
trated in the central column of Fig. 6. This explicit iteration 
method eliminates the large matrix storage requirements of 
steady-state techniques and allows the use of conventional im- 
pedance conditions. As time increases, the iteration process 
directly computes the steady-state variable if/. 

Concluding Remarks 

A transient-frequency domain numerical solution of the po- 
tential acoustic equation has been developed. The potential form 


628 / Vol. 118, OCTOBER 1996 


Transactions of the ASME 




















of the governing equations has been employed to reduce the 
number of dependent variables and their associated storage re- 
quirements. This explicit iteration method represents a signifi- 
cant advance over previous steady-state and transient tech- 
niques. Time is introduced into the steady-state formulation to 
form a hyperbolic equation. A parabolic approximation (in 
time) to this hyperbolic equation is employed in this paper. The 
field is iterated in time from an initial value of 0 to attain the 
steady-state solution. 

The method eliminates the large matrix storage requirements 
of steady state techniques in the frequency domain but still 
allows the use of conventional impedance conditions. Most im- 
portantly, the formulation is explicit under flow conditions. In 
each example provided, the numerical solution quickly and ac- 
curately converges to the exact steady- state solution. 
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